function [AutoCorrYY,StdYY]=VARCovMat(Mat_A,Mat_B,Mat_C,OrderList)
% X_{t+1}=A*X_{t}+B*Shock_{t+1}, Y_{t}=C*X_{t}
% HX=solution.hx;
% M_eta=solution.M_eta;
% GX=solution.gx;
if nargin<4
    OrderList   =   0;
end
Num_X       =   size(Mat_A,1);
Num_S       =   size(Mat_B,2);
Num_Y       =   size(Mat_C,1);
Num_Order   =   length(OrderList);
% Vec_CovXX   =   ( eye(Num_X^2,Num_X^2)-kron(HX,HX) )\VecFun_Mat2Vec(M_eta*M_eta');
% CovXX       =   reshape(Vec_CovXX,[Num_X,Num_X]);


Temp        =   Mat_B*Mat_B';
CovXX       =   Temp;
for tt=1:200
    Temp            =   Mat_A*Temp*Mat_A';
    CovXX           =   CovXX+Temp;
end

CovYY       =   Mat_C*CovXX*Mat_C';
StdYY       =   sqrt(diag(CovYY));
AutoCorrYY  =   zeros(Num_Y,Num_Y,Num_Order);
for ii=1:Num_Order
    AutoCovXX_Temp      =   ( Mat_A^OrderList(ii) )*CovXX;
    AutoCovYY_Temp      =   Mat_C*AutoCovXX_Temp*Mat_C';
    AutoCorrYY(:,:,ii)  =   bsxfun(@rdivide,AutoCovYY_Temp,StdYY*StdYY');
end

if Num_Order==1
    AutoCorrYY      =   squeeze(AutoCorrYY(:,:,1));
end